Sister pair DE analysis, doplots and volcanoplots ctrl -vs- lumb

library(Seurat)
library(dplyr)
library(ggplot2)
library(stringr)
library(tibble)
library(patchwork)
library(plotly)

DE table

First we load the sister pair DE tables and filter for:

• absolute avg_log2FC > 0.5 (~41% increase)

• p_val_adj < 0.01

DE_list <- readRDS("~/spinal_cord_paper/data/Gg_ctrl_lumb_sis_markers.rds")

for (i in seq(DE_list)) {
    DE_list[[i]] <- DE_list[[i]] %>% 
    arrange(desc(avg_log2FC)) %>% 
    filter(abs(avg_log2FC) > 0.5) %>% 
    filter(p_val_adj < 0.01)
}

DE_table <- do.call(rbind, DE_list)
dim(DE_table)

[1] 1510    8

delta pct distribution

par(mfrow = c(2,2))
hist(abs(DE_list[[1]]$delta_pct), breaks = 20)
abline(v = 0.1, lty = "dashed", col = "red")

hist(abs(DE_list[[2]]$delta_pct), breaks = 20)
abline(v = 0.1, lty = "dashed", col = "red")

hist(abs(DE_list[[4]]$delta_pct), breaks = 20)
abline(v = 0.1, lty = "dashed", col = "red")

hist(abs(DE_list[[5]]$delta_pct), breaks = 20)
abline(v = 0.1, lty = "dashed", col = "red")

Now we filter the DE lists for absolute delta percentage > 0.1.

for (i in seq(DE_list)) {
  DE_list[[i]] <- DE_list[[i]] %>% 
  filter(abs(delta_pct) > 0.1)
}

DE_table <- do.call(rbind, DE_list)
dim(DE_table)

[1] 1113    8

Broad clusters

broad_order <- c("progenitors",
      "FP",
      "RP",
      "FP/RP",
      "neurons",
      "OPC",
      "MFOL",
      "pericytes",
      "microglia",
      "blood",
      "vasculature"
      )

Integrated data

Load the integrated control and poly data.

int_path <- "Gg_ctrl_lumb_int_seurat_250723"

my.se <- readRDS(paste0("~/spinal_cord_paper/data/", int_path, ".rds"))
  annot_int <- read.csv(list.files("~/spinal_cord_paper/annotations",
                               pattern = str_remove(int_path, "_seurat_\\d{6}"),
                               full.names = TRUE))
  
  if(length(table(annot_int$number)) != length(table(my.se$seurat_clusters))) {
     stop("Number of clusters must be identical!")
  }
  
  # rename for left join
  annot_int <- annot_int %>% 
    mutate(fine = paste(fine, number, sep = "_")) %>% 
    mutate(number = factor(number, levels = 1:nrow(annot_int))) %>% 
    rename(seurat_clusters = number)
  
  ord_levels <- annot_int$fine[order(match(annot_int$broad, broad_order))]
   
  # add cluster annotation to meta data
  my.se@meta.data <- my.se@meta.data %>% 
    rownames_to_column("rowname") %>% 
    left_join(annot_int, by = "seurat_clusters") %>% 
    mutate(fine = factor(fine, levels = ord_levels)) %>% 
    mutate(seurat_clusters = factor(seurat_clusters, levels = str_extract(ord_levels, "\\d{1,2}$"))) %>% 
    column_to_rownames("rowname")
  
  ctrl_poly_int_combined_labels <- readRDS("~/spinal_cord_paper/annotations/ctrl_lumb_int_combined_labels.rds")
  
  my.se <- AddMetaData(my.se, ctrl_poly_int_combined_labels)
  

DimPlot

DimPlot(
  my.se,
  group.by = "annot_sample",
  reduction = "tsne",
  label = TRUE,
  repel = TRUE
  ) +
  NoLegend()

Warning: ggrepel: 3 unlabeled data points (too many overlaps). Consider increasing max.overlaps

Cluster order

Get the cluster order from the spearman correlation heatmap of the control and poly integrated data. Then we filter for the neuronal clusters only.

corr_heatmap <- readRDS("~/spinal_cord_paper/output/heatmap_spearman_ctrl_lumb.rds")

#heatmap order
htmp_order <- data.frame("label" = corr_heatmap[["gtable"]]$grobs[[4]]$label) %>% 
  mutate(label = str_remove(label, "_int")) %>% 
  mutate(label_ordered = paste(str_sub(label,6 ,-1), str_sub(label, 1, 4), sep = "_"))

my.se@meta.data <- my.se@meta.data %>%
  mutate(annot_sample = factor(annot_sample, levels = htmp_order$label_ordered))

Idents(my.se) <- "annot_sample"

# filter for the neuronal clusters
my.se <- subset(my.se, idents = htmp_order$label_ordered[grepl("neurons|MN", htmp_order$label_ordered)])

DimPlot(
  my.se,
  group.by = "annot_sample",
  reduction = "tsne",
  label = TRUE,
  repel = TRUE
  ) +
  NoLegend()

my.se@active.assay <- "RNA"

Individual dot plots

# select top50 by log2FC 
for (i in seq(DE_list)) {
    DE_list[[i]] <- DE_list[[i]] %>%
    slice_max(order_by = abs(avg_log2FC), n = 50) %>% 
    arrange(desc(avg_log2FC))
}

p1 <- modplots::mDotPlot2(my.se,
                    group.by = "annot_sample", 
                    assay = "RNA",
                      # reverse order of DE genes so number one is on top
                    features = rev(DE_list[[1]]$Gene.stable.ID),
                    gnames = modplots::gnames,
          cols = c("lightgrey", "black")) +
    theme(axis.text.x = element_text(angle = 90, hjust=1, vjust=0.5)) +
    coord_flip() +
    xlab(names(DE_list)[1])

p2 <- modplots::mDotPlot2(my.se,
                    group.by = "annot_sample",  
                    assay = "RNA",
                      # reverse order of DE genes so number one is on top
                    features = rev(DE_list[[2]]$Gene.stable.ID),
                    gnames = modplots::gnames,
          cols = c("lightgrey", "black")) +
    theme(axis.text.x = element_text(angle = 90, hjust=1, vjust=0.5)) +
    coord_flip() +
    xlab(names(DE_list)[2])

p3 <- modplots::mDotPlot2(my.se,
                    group.by = "annot_sample",  
                    assay = "RNA",
                      # reverse order of DE genes so number one is on top
                    features = rev(DE_list[[3]]$Gene.stable.ID),
                    gnames = modplots::gnames,
          cols = c("lightgrey", "black")) +
    theme(axis.text.x = element_text(angle = 90, hjust=1, vjust=0.5)) +
    coord_flip() +
    xlab(names(DE_list)[3])

p4 <- modplots::mDotPlot2(my.se,
                    group.by = "annot_sample",  
                    assay = "RNA",
                      # reverse order of DE genes so number one is on top
                    features = rev(DE_list[[4]]$Gene.stable.ID),
                    gnames = modplots::gnames,
          cols = c("lightgrey", "black")) +
    theme(axis.text.x = element_text(angle = 90, hjust=1, vjust=0.5)) +
    coord_flip() +
    xlab(names(DE_list)[4])

p5 <- modplots::mDotPlot2(my.se,
                    group.by = "annot_sample",  
                    assay = "RNA",
                      # reverse order of DE genes so number one is on top
                    features = rev(DE_list[[5]]$Gene.stable.ID),
                    gnames = modplots::gnames,
          cols = c("lightgrey", "black")) +
    theme(axis.text.x = element_text(angle = 90, hjust=1, vjust=0.5)) +
    coord_flip() +
    xlab(names(DE_list)[5])

pdf("~/spinal_cord_paper/figures/ctrl_lumb_dotplot_individual.pdf", height = 13, width = 20)
(p1 + p2 + p3 + p4 + p5) + plot_layout(guides = "collect", nrow = 1)
dev.off()

null device 
          1 

Volcanoplots

p.adj <- 0.01
l2fc <- 0.5

# select top50 by log2FC 
for (i in seq(DE_list)) {
    DE_list[[i]] <- DE_list[[i]] %>% 
    mutate(delta_pct_sign = case_when(
      delta_pct < 0 ~ "-",
      delta_pct > 0 ~ "+",
      delta_pct == 0 ~ "0"
    ))
}
 

toplot <- do.call(rbind, DE_list) %>% 
  rownames_to_column("contrast") %>% 
  mutate(contrast = str_remove(contrast, "\\.\\d{1,2}")) %>% 
  mutate(contrast = str_replace_all(contrast, " ", "_")) %>% 
  filter(!grepl("^HOX", Gene.name)) # remove hox genes

volplot <- ggplot(data = toplot,
       aes(x = avg_log2FC,
           y = -log10(p_val_adj),
           label = Gene.name,
           color = delta_pct_sign,
           size = abs(delta_pct)
       )) +
  geom_point(shape = 21) +
  geom_hline(yintercept = -log10(p.adj), linetype = "dashed") +
  geom_vline(xintercept = c(-l2fc,l2fc), linetype = "dashed") +
  scale_color_manual(values = c("#419c73", "black")) +
  scale_size_continuous(range = c(0.5, 4)) +
  facet_wrap("contrast", ncol = 5, scales = "free") +
  ylab("-log10(padj)") +
  theme_bw()

ggplotly(volplot)

NA

pdf("~/spinal_cord_paper/figures/Fig_4_volcanoplots.pdf", width = 15, height = 15)
(volplot +
  ggrepel::geom_text_repel(size = 3, color = "black"))

Warning: ggrepel: 5 unlabeled data points (too many overlaps). Consider increasing max.overlaps
Warning: ggrepel: 11 unlabeled data points (too many overlaps). Consider increasing max.overlaps
Warning: ggrepel: 7 unlabeled data points (too many overlaps). Consider increasing max.overlaps
Warning: ggrepel: 6 unlabeled data points (too many overlaps). Consider increasing max.overlaps
Warning: ggrepel: 6 unlabeled data points (too many overlaps). Consider increasing max.overlaps
Warning: ggrepel: 5 unlabeled data points (too many overlaps). Consider increasing max.overlaps
Warning: ggrepel: 29 unlabeled data points (too many overlaps). Consider increasing max.overlaps

# Date and time of Rendering
Sys.time()

sessionInfo()